##############################################################################

library(CHAT)
library(argparser)

##############################################################################

argp <- arg_parser("Different ways of handling SNVs not adjusted by CHAT and SNVs in A1 or A2 lineage")
argp <- add_argument(argp, "--mCRC_dir", help="path to the mCRC dir")
argp <- add_argument(argp, "--Sample", help="Sample")
argp <- add_argument(argp, "--mut_file", help="Mut file , merge all samples from one people")
argp <- add_argument(argp, "--seg_dir", help="titan cnv path")
argp <- add_argument(argp, "--purity_file", help="purity_file")
argp <- add_argument(argp, "--out_path", help="Output path")

argv <- parse_args(argp)

if(1!=1){
  mCRC_dir <- "/public/home/xxf2019/tools/mCRCs"
  mut_file <- "~/20220915_gastric_multiple/dna_combine/titan/mut/S1_Mut.tsv"
  seg_dir <- "~/20220915_gastric_multiple/dna_combine/titan/FinalModel"
  out_path <- "~/20220915_gastric_multiple/dna_combine/titan/chat"
  purity_file <- "~/20220915_gastric_multiple/dna_combine/titan/Purity_titan.final.tsv"
  Sample <- "S1"
}

mCRC_dir <- argv$mCRC_dir
mut_file <- argv$mut_file
seg_dir <- argv$seg_dir
Sample <- argv$Sample
out_path <- argv$out_path
purity_file <- argv$purity_file

##############################################################################

dir1 <- paste0(mCRC_dir,"/data_processing/")
source(paste0(dir1,"getSampleCCF2.R"))
source(paste0(dir1,"TitanCNA_04_TitanCNA2seg.R"))

##############################################################################

#chrs <- paste0("chr",c(1:22))
chrs <- c(1:22)

count1 <- read.delim(mut_file,as.is=TRUE)
count1 <- count1[count1$Chr %in% chrs,]
count1$Pos <- count1$Start_Position

Tumor <- Sample

#############################
## 读取CNV的文件
## titan最终输出的cnv文件为后缀是.ichor.seg.txt
# cnv_path <- paste0(seg_dir,"/", Tumor ,"_titanCNA_result")
# cnv_files <- list.files(cnv_path)
cnv_files <- list.files(seg_dir)
seg_file <- grep( "seg" , grep(paste0("^" , Sample , "_") , cnv_files , value=T) , value = T)

dat_seg <- read.table(paste0(seg_dir,"/",seg_file) , header = T)
dat_purity <- read.table( purity_file , header = T )
dat_params <- subset( dat_purity , Sample == Tumor )
colnames(dat_params) <- c("Tumor" , "purity" , "ploidy")

cnaseg <- data.frame(chrom=dat_seg$Chromosome ,
                  loc.start=dat_seg$Start ,
                  loc.end=dat_seg$End ,
                  num.mark=dat_seg$Length.snp. ,
                  seg.mean=dat_seg$Corrected_logR,
      copynumber=dat_seg$Corrected_Copy_Number,
      minor_cn=dat_seg$Corrected_MinorCN,
      major_cn=dat_seg$Corrected_MajorCN,
      allelicratio=dat_seg$Corrected_Ratio,
      LOHcall=dat_seg$Corrected_Call,
      cellularprevalence=dat_seg$Cellular_Prevalence,
      ploidy=dat_params$ploidy,
      normalproportion=(1-dat_params$purity)
    )

cna1 <- cnaseg[,c(1,2,3,5,4,9,4,11,7,6)]
cna1[,8] <- cna1[,8]*(1-cnaseg$normalproportion[1])
cna1[,1] <- as.character(cna1[,1])
cna1$cellularprevalence[is.na(cna1$cellularprevalence)] <- max(cna1$cellularprevalence,na.rm=TRUE)
cna1 <- cna1[cna1$chrom %in% chrs,]

############################
cover1 <- count1$cover
freq1 <- count1$freq

ref1 <- round(cover1 * (1-freq1))
var1 <- round(cover1 * freq1)
snv1 <- data.frame(Chr=count1$Chr,Pos=count1$Start_Position,
	     Tumor_cover=cover1,Tumor_freq=freq1,
	     stringsAsFactors=FALSE)
############################

############################
## 计算CCF
ccf1 <- getSampleCCF2(snv1,cna1)
## 无CNV的样本，默认处于二倍体区域，CCF=2*VAF
ccf1$CCF[is.na(ccf1$CCF)] <- 2*ccf1$Tumor_freq[is.na(ccf1$CCF)]
ccf1$CCF[ccf1$Tumor_freq==0] <- 0

ccf1$LOH <- rep(FALSE,nrow(ccf1))
ccf1$LOH[ccf1$n_minor == 0 & ccf1$sAGP == max(ccf1$sAGP)] <- TRUE

purity1 <- dat_params$purity
ccf1a <- ccf1$CCF/purity1

count1 <- cbind(count1,round(2*freq1,3),round(2*freq1/purity1,3),FALSE)
colnames(count1)[ncol(count1)-(2:0)] <- c("CCF","CCF_adj","LOH")
idx1 <- match(paste(ccf1$Chr,ccf1$Pos,sep="_"),paste(count1$Chr,count1$Pos,sep="_"))
count1[idx1,ncol(count1)-2] <- round(ccf1$CCF,3)
count1[idx1,ncol(count1)-1] <- round(ccf1a,3)
count1[idx1,ncol(count1)] <- ccf1$LOH

##### 单个样本的输出
idx1 <- match(paste(ccf1$Chr,ccf1$Pos,sep="_"),paste(count1$Chr,count1$Start_Position,sep="_"))
count1[idx1,"minor_cn"] <- ccf1$n_minor
count1[idx1,"total_cn"] <- ccf1$n_total
##### 修改对应的列名
colnames(count1) <- c("Chr" , "Start_Position" , "End_Position" , "REF" , "ALT" , "Hugo_Symbol" , "Variant_Classification" ,
    "Depth" , "VAF" , "ID" , "Pos" , "CCF" , "CCF_adj" , "LOH" , "minor_cn" , "total_cn"
  )
count1 <- count1[,c("Chr" , "Start_Position" , "End_Position" , 
  "REF" , "ALT" , "Hugo_Symbol" , "Variant_Classification" ,
  "Depth" , "VAF" , "CCF_adj" , "minor_cn" , "total_cn")]

##truncate at 1
count1$CCF_adj <- ifelse(count1$CCF_adj > 1 , 1 , count1$CCF_adj)
## 去除无拷贝覆盖的区域
## count1 <- count1[!(is.na(count1$minor_cn)),]

##### 输出
out_file <- paste0( out_path , "/" , Tumor  , "_CHAT.txt")
write.table(count1 , out_file, quote=FALSE,sep="\t",row.names=F)
